function [ResEuler,xGrid] = EulerEqError_Proj_mp(model,setup,orderApp)

% Some indices 
ny      = model.ny;        %Number of output variables
nx      = model.nx;        %Number of state variables 
 
%Unfold params
BETTA= model.params.BETTA;
B= model.params.B;
CHI= model.params.CHI;
CHI0= model.params.CHI0;
THETA= model.params.THETA;
DELTA= model.params.DELTA;
ALFA= model.params.ALFA;
PHI= model.params.PHI;
PHIzero= model.params.PHIzero;
ZI= model.params.ZI;
KAPAw= model.params.KAPAw;
ETA= model.params.ETA;
RHOR= model.params.RHOR;
PHIpai= model.params.PHIpai;
PHIpai_1= model.params.PHIpai_1;
PHIy= model.params.PHIy;
PHIy_1= model.params.PHIy_1;
PHIc= model.params.PHIc;
PHIc_1= model.params.PHIc_1;
PHIl= model.params.PHIl;
NU= model.params.NU;
U0= model.params.U0;
U0d= model.params.U0d;
OMEGAz= model.params.OMEGAz;
OMEGAd= model.params.OMEGAd;
OMEGAn= model.params.OMEGAn;
OMEGAp= model.params.OMEGAp;
OMEGAa= model.params.OMEGAa;
RHOz= model.params.RHOz;
RHOd= model.params.RHOd;
RHOn= model.params.RHOn;
RHOp= model.params.RHOp;
RHOa= model.params.RHOa;
Kss= model.params.Kss;
OUTPUTss= model.params.OUTPUTss;
PAIss= model.params.PAIss;
MUZss= model.params.MUZss;
Css= model.params.Css;
AA= model.params.AA;
lss= model.params.lss;
Rss= model.params.Rss;
Wss= model.params.Wss;
 
% The raw grid for the states 
if setup.xPointsDim == 0 
    xGrid = transpose(setup.xSim); 
else 
    dx       = (setup.upperx-setup.lowerx)/(setup.xPointsDim-1); 
    xGridtmp = zeros(nx,setup.xPointsDim); 
    for j=1:setup.xPointsDim 
        xGridtmp(:,j) = setup.lowerx(:,1)+(j-1)*dx(:,1); 
    end 
    tmp = cell(1,nx); 
    for i=1:nx 
        tmp(i) = mat2cell(xGridtmp(i,:),1,setup.xPointsDim); 
    end 
    xGrid = cartprod(tmp); 
end 
 
% The grid for numerical itegration 
tmp_e        = cell(1,model.ne); 
tmp_w        = cell(1,model.ne); 
for i=1:model.ne 
    tmp_e(i) = mat2cell(setup.GH_e,length(setup.GH_e),1); 
    tmp_w(i) = mat2cell(setup.GH_w,length(setup.GH_w),1); 
end 
eGrid = cartprod(tmp_e); 
wGrid = prod(cartprod(tmp_w),2); 
 
%% Building the residuals 
ResEuler = zeros(size(xGrid,1),ny+nx); 
 
%start of multiprocessing
numCPUs = feature('numcores');
parpool('local',numCPUs);
parfor (k=1:size(xGrid,1),numCPUs)    
    %if mod(k,1000) == 0
    %    disp(['Points done so far = ',num2str(k)])
    %end
    x_cu   = transpose(xGrid(k,:)); 
    % x_cu = x_t - xss
    % Getting the states at time t 
    c_ba1 = model.hSS(1,1) + x_cu(1,1);
    muz_cu = model.hSS(2,1) + x_cu(2,1);
    d_cu = model.hSS(3,1) + x_cu(3,1);
    n_cu = model.hSS(4,1) + x_cu(4,1);
    paistar_cu = model.hSS(5,1) + x_cu(5,1);
    a_cu = model.hSS(6,1) + x_cu(6,1);
 
    % Getting the controls at time t and part of the states at time t+1 
    if orderApp == 1
        yVecA_cu = model.g0  + model.gx*x_cu;
        % Bond yields to bond prices
        scalar   = -(1:40)';
        yVecB_cu = model.byx0 + model.by*x_cu;
        yVec_cu  = [yVecA_cu; scalar.*yVecB_cu];
        
        xVec_cup = model.hSS + model.h0 + model.hx*x_cu;
    elseif orderApp == 2
        XX_cu = kron(x_cu,x_cu);
        yVecA_cu = model.g0  + model.gx*x_cu + model.Gxx*XX_cu;
        % Bond yields to bond prices
        scalar   = -(1:40)';
        yVecB_cu = model.by0 + model.byx*x_cu + model.Byxx*XX_cu;
        yVec_cu  = [yVecA_cu; scalar.*yVecB_cu];
        
        xVec_cup = model.hSS + model.h0 + model.hx*x_cu + model.Hxx*XX_cu;
    elseif orderApp == 3
        XX_cu  = kron(x_cu,x_cu);
        XXX_cu = kron(XX_cu,x_cu);
        yVecA_cu = model.g0  + model.gx*x_cu + model.Gxx*XX_cu + model.Gxxx*XXX_cu;
        % Bond yields to bond prices
        scalar   = -(1:40)';
        yVecB_cu = model.by0 + model.byx*x_cu + model.Byxx*XX_cu + model.Byxxx*XXX_cu;
        yVec_cu  = [yVecA_cu; scalar.*yVecB_cu];
        
        xVec_cup = model.hSS + model.h0 + model.hx*x_cu + model.Hxx*XX_cu + model.Hxxx*XXX_cu;

    end

    % The controls 
    c_cu = yVec_cu(1,1);
    r_cu = yVec_cu(2,1);
    pai_cu = yVec_cu(3,1);
    w_cu = yVec_cu(4,1);
    l_cu = yVec_cu(5,1);
    output_cu = yVec_cu(6,1);
    mc_cu = yVec_cu(7,1);
    evf_cu = yVec_cu(8,1);
    vf_cu = yVec_cu(9,1);
    p1_cu = yVec_cu(10,1);
    p2_cu = yVec_cu(11,1);
    p3_cu = yVec_cu(12,1);
    p4_cu = yVec_cu(13,1);
    p5_cu = yVec_cu(14,1);
    p6_cu = yVec_cu(15,1);
    p7_cu = yVec_cu(16,1);
    p8_cu = yVec_cu(17,1);
    p9_cu = yVec_cu(18,1);
    p10_cu = yVec_cu(19,1);
    p11_cu = yVec_cu(20,1);
    p12_cu = yVec_cu(21,1);
    p13_cu = yVec_cu(22,1);
    p14_cu = yVec_cu(23,1);
    p15_cu = yVec_cu(24,1);
    p16_cu = yVec_cu(25,1);
    p17_cu = yVec_cu(26,1);
    p18_cu = yVec_cu(27,1);
    p19_cu = yVec_cu(28,1);
    p20_cu = yVec_cu(29,1);
    p21_cu = yVec_cu(30,1);
    p22_cu = yVec_cu(31,1);
    p23_cu = yVec_cu(32,1);
    p24_cu = yVec_cu(33,1);
    p25_cu = yVec_cu(34,1);
    p26_cu = yVec_cu(35,1);
    p27_cu = yVec_cu(36,1);
    p28_cu = yVec_cu(37,1);
    p29_cu = yVec_cu(38,1);
    p30_cu = yVec_cu(39,1);
    p31_cu = yVec_cu(40,1);
    p32_cu = yVec_cu(41,1);
    p33_cu = yVec_cu(42,1);
    p34_cu = yVec_cu(43,1);
    p35_cu = yVec_cu(44,1);
    p36_cu = yVec_cu(45,1);
    p37_cu = yVec_cu(46,1);
    p38_cu = yVec_cu(47,1);
    p39_cu = yVec_cu(48,1);
    p40_cu = yVec_cu(49,1);
 
    Euler = zeros(ny+nx,1);
    for i=1:length(wGrid)
        x_cup = xVec_cup-model.hSS; 
        % Adding shocks to the states
        x_cup = x_cup + model.eta*transpose(eGrid(i,:));
        % Getting controls at time t+1 
        if orderApp == 1
            yVecA_cup = model.g0  + model.gx*x_cup;
            % Bond yields to bond prices
            scalar    = -(1:40)';
            yVecB_cup = model.by0 + model.byx*x_cup;
            yVec_cup  = [yVecA_cup; scalar.*yVecB_cup];
        elseif orderApp == 2
            XX_cup    = kron(x_cup,x_cup);
            yVecA_cup = model.g0  + model.gx*x_cup + model.Gxx*XX_cup;
            % Bond yields to bond prices
            scalar    = -(1:40)';
            yVecB_cup = model.by0 + model.byx*x_cup + model.Byxx*XX_cup;
            yVec_cup  = [yVecA_cup; scalar.*yVecB_cup];
            
        elseif orderApp == 3
            XX_cup    = kron(x_cup,x_cup);
            XXX_cup   = kron(XX_cup,x_cup);
            yVecA_cup = model.g0  + model.gx*x_cup + model.Gxx*XX_cup + model.Gxxx*XXX_cup;
            % Bond yields to bond prices
            scalar    = -(1:40)';
            yVecB_cup = model.by0 + model.byx*x_cup + model.Byxx*XX_cup + model.Byxxx*XXX_cup;
            yVec_cup  = [yVecA_cup; scalar.*yVecB_cup];
        end
        
        c_cup = yVec_cup(1,1);
        r_cup = yVec_cup(2,1);
        pai_cup = yVec_cup(3,1);
        w_cup = yVec_cup(4,1);
        l_cup = yVec_cup(5,1);
        output_cup = yVec_cup(6,1);
        mc_cup = yVec_cup(7,1);
        evf_cup = yVec_cup(8,1);
        vf_cup = yVec_cup(9,1);
        p1_cup = yVec_cup(10,1);
        p2_cup = yVec_cup(11,1);
        p3_cup = yVec_cup(12,1);
        p4_cup = yVec_cup(13,1);
        p5_cup = yVec_cup(14,1);
        p6_cup = yVec_cup(15,1);
        p7_cup = yVec_cup(16,1);
        p8_cup = yVec_cup(17,1);
        p9_cup = yVec_cup(18,1);
        p10_cup = yVec_cup(19,1);
        p11_cup = yVec_cup(20,1);
        p12_cup = yVec_cup(21,1);
        p13_cup = yVec_cup(22,1);
        p14_cup = yVec_cup(23,1);
        p15_cup = yVec_cup(24,1);
        p16_cup = yVec_cup(25,1);
        p17_cup = yVec_cup(26,1);
        p18_cup = yVec_cup(27,1);
        p19_cup = yVec_cup(28,1);
        p20_cup = yVec_cup(29,1);
        p21_cup = yVec_cup(30,1);
        p22_cup = yVec_cup(31,1);
        p23_cup = yVec_cup(32,1);
        p24_cup = yVec_cup(33,1);
        p25_cup = yVec_cup(34,1);
        p26_cup = yVec_cup(35,1);
        p27_cup = yVec_cup(36,1);
        p28_cup = yVec_cup(37,1);
        p29_cup = yVec_cup(38,1);
        p30_cup = yVec_cup(39,1);
        p31_cup = yVec_cup(40,1);
        p32_cup = yVec_cup(41,1);
        p33_cup = yVec_cup(42,1);
        p34_cup = yVec_cup(43,1);
        p35_cup = yVec_cup(44,1);
        p36_cup = yVec_cup(45,1);
        p37_cup = yVec_cup(46,1);
        p38_cup = yVec_cup(47,1);
        p39_cup = yVec_cup(48,1);
        p40_cup = yVec_cup(49,1);
 
        % Getting the states at time t+1 
        c_ba1p = model.hSS(1,1) + x_cup(1,1);
        muz_cup = model.hSS(2,1) + x_cup(2,1);
        d_cup = model.hSS(3,1) + x_cup(3,1);
        n_cup = model.hSS(4,1) + x_cup(4,1);
        paistar_cup = model.hSS(5,1) + x_cup(5,1);
        a_cup = model.hSS(6,1) + x_cup(6,1);
        %% Model Specific 
        f     = zeros(ny+nx,1);
% START DISPLAYING f
f(1,1)=exp(-vf_cu)*(exp(d_cu)*(((exp(c_cu) - B*exp(-muz_cu)*exp(c_ba1))/Css^CHI0)^(1 - CHI)/(CHI - 1) - U0d + (PHIzero*exp(n_cu)*(1 - exp(l_cu))^(1 - 1/PHI))/(1/PHI - 1)) - U0 + (AA*BETTA)/exp(evf_cu)^(1/(ALFA - 1))) - 1;
f(2,1)=exp(-w_cu)*(KAPAw*Wss - (Css^CHI0*PHIzero*exp(n_cu)*((exp(c_cu) - B*exp(-muz_cu)*exp(c_ba1))/Css^CHI0)^CHI*(KAPAw - 1))/(1 - exp(l_cu))^(1/PHI)) - 1;
f(3,1)=(BETTA*exp(-d_cu)*exp(-pai_cup)*exp(d_cup)*exp(r_cu)*exp(muz_cup)^(CHI*(CHI0 - 1) - CHI0)*(exp(c_cu) - B*exp(-muz_cu)*exp(c_ba1))^CHI*((AA*exp(-vf_cup))/(exp(evf_cu)^(1/(ALFA - 1))*exp(muz_cup)^((CHI - 1)*(CHI0 - 1))))^ALFA)/(exp(c_cup) - B*exp(-muz_cup)*exp(c_cu))^CHI - 1;
f(4,1)=- (exp(-a_cu)*exp(-mc_cu)*exp(w_cu)*exp(l_cu)^THETA)/(Kss^THETA*(THETA - 1)) - 1;
f(5,1)=1 - (exp(-mc_cu)*(ETA + (ZI*exp(pai_cu)*(exp(pai_cu)/PAIss^NU - 1))/PAIss^NU - (BETTA*ZI*exp(-d_cu)*exp(-output_cu)*exp(d_cup)*exp(muz_cup)*exp(output_cup)*exp(pai_cup)*exp(muz_cup)^(CHI*(CHI0 - 1) - CHI0)*(exp(c_cu) - B*exp(-muz_cu)*exp(c_ba1))^CHI*(exp(pai_cup)/PAIss^NU - 1)*((AA*exp(-vf_cup))/(exp(evf_cu)^(1/(ALFA - 1))*exp(muz_cup)^((CHI - 1)*(CHI0 - 1))))^ALFA)/(PAIss^NU*(exp(c_cup) - B*exp(-muz_cup)*exp(c_cu))^CHI) - 1))/ETA;
f(6,1)=Rss*exp(PHIl*log(exp(l_cu)/lss) + PHIc*(log(exp(muz_cu)) - log(MUZss) - exp(c_ba1) + exp(c_cu)) + PHIc_1*(log(exp(muz_cup)) - log(MUZss) - exp(c_ba1p) + exp(c_cup)) - PHIpai*(log(exp(paistar_cu)) - log(exp(pai_cu)) + log(PAIss)) - PHIpai_1*(log(exp(paistar_cup)) - log(exp(pai_cup)) + log(PAIss)) + PHIy*log(exp(output_cu)/OUTPUTss) + PHIy_1*log(exp(output_cup)/OUTPUTss)) - exp(r_cu);
f(7,1)=(exp(output_cu)*((ZI*(exp(pai_cu)/PAIss^NU - 1)^2)/2 - 1))/(exp(c_cu) + DELTA*Kss) + 1;
f(8,1)=Kss^THETA*exp(-output_cu)*exp(a_cu)*exp(l_cu)^(1 - THETA) - 1;
f(9,1)=exp(-evf_cu)*((exp(vf_cup)*exp(muz_cup)^((CHI - 1)*(CHI0 - 1)))/AA)^(1 - ALFA) - 1;
f(10,1)=exp(-p1_cu)*exp(-r_cu) - 1;
f(11,1)=(BETTA*exp(-d_cu)*exp(-p2_cu)*exp(-pai_cup)*exp(d_cup)*exp(p1_cup)*exp(muz_cup)^(CHI*(CHI0 - 1) - CHI0)*(exp(c_cu) - B*exp(-muz_cu)*exp(c_ba1))^CHI*((AA*exp(-vf_cup))/(exp(evf_cu)^(1/(ALFA - 1))*exp(muz_cup)^((CHI - 1)*(CHI0 - 1))))^ALFA)/(exp(c_cup) - B*exp(-muz_cup)*exp(c_cu))^CHI - 1;
f(12,1)=(BETTA*exp(-d_cu)*exp(-p3_cu)*exp(-pai_cup)*exp(d_cup)*exp(p2_cup)*exp(muz_cup)^(CHI*(CHI0 - 1) - CHI0)*(exp(c_cu) - B*exp(-muz_cu)*exp(c_ba1))^CHI*((AA*exp(-vf_cup))/(exp(evf_cu)^(1/(ALFA - 1))*exp(muz_cup)^((CHI - 1)*(CHI0 - 1))))^ALFA)/(exp(c_cup) - B*exp(-muz_cup)*exp(c_cu))^CHI - 1;
f(13,1)=(BETTA*exp(-d_cu)*exp(-p4_cu)*exp(-pai_cup)*exp(d_cup)*exp(p3_cup)*exp(muz_cup)^(CHI*(CHI0 - 1) - CHI0)*(exp(c_cu) - B*exp(-muz_cu)*exp(c_ba1))^CHI*((AA*exp(-vf_cup))/(exp(evf_cu)^(1/(ALFA - 1))*exp(muz_cup)^((CHI - 1)*(CHI0 - 1))))^ALFA)/(exp(c_cup) - B*exp(-muz_cup)*exp(c_cu))^CHI - 1;
f(14,1)=(BETTA*exp(-d_cu)*exp(-p5_cu)*exp(-pai_cup)*exp(d_cup)*exp(p4_cup)*exp(muz_cup)^(CHI*(CHI0 - 1) - CHI0)*(exp(c_cu) - B*exp(-muz_cu)*exp(c_ba1))^CHI*((AA*exp(-vf_cup))/(exp(evf_cu)^(1/(ALFA - 1))*exp(muz_cup)^((CHI - 1)*(CHI0 - 1))))^ALFA)/(exp(c_cup) - B*exp(-muz_cup)*exp(c_cu))^CHI - 1;
f(15,1)=(BETTA*exp(-d_cu)*exp(-p6_cu)*exp(-pai_cup)*exp(d_cup)*exp(p5_cup)*exp(muz_cup)^(CHI*(CHI0 - 1) - CHI0)*(exp(c_cu) - B*exp(-muz_cu)*exp(c_ba1))^CHI*((AA*exp(-vf_cup))/(exp(evf_cu)^(1/(ALFA - 1))*exp(muz_cup)^((CHI - 1)*(CHI0 - 1))))^ALFA)/(exp(c_cup) - B*exp(-muz_cup)*exp(c_cu))^CHI - 1;
f(16,1)=(BETTA*exp(-d_cu)*exp(-p7_cu)*exp(-pai_cup)*exp(d_cup)*exp(p6_cup)*exp(muz_cup)^(CHI*(CHI0 - 1) - CHI0)*(exp(c_cu) - B*exp(-muz_cu)*exp(c_ba1))^CHI*((AA*exp(-vf_cup))/(exp(evf_cu)^(1/(ALFA - 1))*exp(muz_cup)^((CHI - 1)*(CHI0 - 1))))^ALFA)/(exp(c_cup) - B*exp(-muz_cup)*exp(c_cu))^CHI - 1;
f(17,1)=(BETTA*exp(-d_cu)*exp(-p8_cu)*exp(-pai_cup)*exp(d_cup)*exp(p7_cup)*exp(muz_cup)^(CHI*(CHI0 - 1) - CHI0)*(exp(c_cu) - B*exp(-muz_cu)*exp(c_ba1))^CHI*((AA*exp(-vf_cup))/(exp(evf_cu)^(1/(ALFA - 1))*exp(muz_cup)^((CHI - 1)*(CHI0 - 1))))^ALFA)/(exp(c_cup) - B*exp(-muz_cup)*exp(c_cu))^CHI - 1;
f(18,1)=(BETTA*exp(-d_cu)*exp(-p9_cu)*exp(-pai_cup)*exp(d_cup)*exp(p8_cup)*exp(muz_cup)^(CHI*(CHI0 - 1) - CHI0)*(exp(c_cu) - B*exp(-muz_cu)*exp(c_ba1))^CHI*((AA*exp(-vf_cup))/(exp(evf_cu)^(1/(ALFA - 1))*exp(muz_cup)^((CHI - 1)*(CHI0 - 1))))^ALFA)/(exp(c_cup) - B*exp(-muz_cup)*exp(c_cu))^CHI - 1;
f(19,1)=(BETTA*exp(-d_cu)*exp(-p10_cu)*exp(-pai_cup)*exp(d_cup)*exp(p9_cup)*exp(muz_cup)^(CHI*(CHI0 - 1) - CHI0)*(exp(c_cu) - B*exp(-muz_cu)*exp(c_ba1))^CHI*((AA*exp(-vf_cup))/(exp(evf_cu)^(1/(ALFA - 1))*exp(muz_cup)^((CHI - 1)*(CHI0 - 1))))^ALFA)/(exp(c_cup) - B*exp(-muz_cup)*exp(c_cu))^CHI - 1;
f(20,1)=(BETTA*exp(-d_cu)*exp(-p11_cu)*exp(-pai_cup)*exp(d_cup)*exp(p10_cup)*exp(muz_cup)^(CHI*(CHI0 - 1) - CHI0)*(exp(c_cu) - B*exp(-muz_cu)*exp(c_ba1))^CHI*((AA*exp(-vf_cup))/(exp(evf_cu)^(1/(ALFA - 1))*exp(muz_cup)^((CHI - 1)*(CHI0 - 1))))^ALFA)/(exp(c_cup) - B*exp(-muz_cup)*exp(c_cu))^CHI - 1;
f(21,1)=(BETTA*exp(-d_cu)*exp(-p12_cu)*exp(-pai_cup)*exp(d_cup)*exp(p11_cup)*exp(muz_cup)^(CHI*(CHI0 - 1) - CHI0)*(exp(c_cu) - B*exp(-muz_cu)*exp(c_ba1))^CHI*((AA*exp(-vf_cup))/(exp(evf_cu)^(1/(ALFA - 1))*exp(muz_cup)^((CHI - 1)*(CHI0 - 1))))^ALFA)/(exp(c_cup) - B*exp(-muz_cup)*exp(c_cu))^CHI - 1;
f(22,1)=(BETTA*exp(-d_cu)*exp(-p13_cu)*exp(-pai_cup)*exp(d_cup)*exp(p12_cup)*exp(muz_cup)^(CHI*(CHI0 - 1) - CHI0)*(exp(c_cu) - B*exp(-muz_cu)*exp(c_ba1))^CHI*((AA*exp(-vf_cup))/(exp(evf_cu)^(1/(ALFA - 1))*exp(muz_cup)^((CHI - 1)*(CHI0 - 1))))^ALFA)/(exp(c_cup) - B*exp(-muz_cup)*exp(c_cu))^CHI - 1;
f(23,1)=(BETTA*exp(-d_cu)*exp(-p14_cu)*exp(-pai_cup)*exp(d_cup)*exp(p13_cup)*exp(muz_cup)^(CHI*(CHI0 - 1) - CHI0)*(exp(c_cu) - B*exp(-muz_cu)*exp(c_ba1))^CHI*((AA*exp(-vf_cup))/(exp(evf_cu)^(1/(ALFA - 1))*exp(muz_cup)^((CHI - 1)*(CHI0 - 1))))^ALFA)/(exp(c_cup) - B*exp(-muz_cup)*exp(c_cu))^CHI - 1;
f(24,1)=(BETTA*exp(-d_cu)*exp(-p15_cu)*exp(-pai_cup)*exp(d_cup)*exp(p14_cup)*exp(muz_cup)^(CHI*(CHI0 - 1) - CHI0)*(exp(c_cu) - B*exp(-muz_cu)*exp(c_ba1))^CHI*((AA*exp(-vf_cup))/(exp(evf_cu)^(1/(ALFA - 1))*exp(muz_cup)^((CHI - 1)*(CHI0 - 1))))^ALFA)/(exp(c_cup) - B*exp(-muz_cup)*exp(c_cu))^CHI - 1;
f(25,1)=(BETTA*exp(-d_cu)*exp(-p16_cu)*exp(-pai_cup)*exp(d_cup)*exp(p15_cup)*exp(muz_cup)^(CHI*(CHI0 - 1) - CHI0)*(exp(c_cu) - B*exp(-muz_cu)*exp(c_ba1))^CHI*((AA*exp(-vf_cup))/(exp(evf_cu)^(1/(ALFA - 1))*exp(muz_cup)^((CHI - 1)*(CHI0 - 1))))^ALFA)/(exp(c_cup) - B*exp(-muz_cup)*exp(c_cu))^CHI - 1;
f(26,1)=(BETTA*exp(-d_cu)*exp(-p17_cu)*exp(-pai_cup)*exp(d_cup)*exp(p16_cup)*exp(muz_cup)^(CHI*(CHI0 - 1) - CHI0)*(exp(c_cu) - B*exp(-muz_cu)*exp(c_ba1))^CHI*((AA*exp(-vf_cup))/(exp(evf_cu)^(1/(ALFA - 1))*exp(muz_cup)^((CHI - 1)*(CHI0 - 1))))^ALFA)/(exp(c_cup) - B*exp(-muz_cup)*exp(c_cu))^CHI - 1;
f(27,1)=(BETTA*exp(-d_cu)*exp(-p18_cu)*exp(-pai_cup)*exp(d_cup)*exp(p17_cup)*exp(muz_cup)^(CHI*(CHI0 - 1) - CHI0)*(exp(c_cu) - B*exp(-muz_cu)*exp(c_ba1))^CHI*((AA*exp(-vf_cup))/(exp(evf_cu)^(1/(ALFA - 1))*exp(muz_cup)^((CHI - 1)*(CHI0 - 1))))^ALFA)/(exp(c_cup) - B*exp(-muz_cup)*exp(c_cu))^CHI - 1;
f(28,1)=(BETTA*exp(-d_cu)*exp(-p19_cu)*exp(-pai_cup)*exp(d_cup)*exp(p18_cup)*exp(muz_cup)^(CHI*(CHI0 - 1) - CHI0)*(exp(c_cu) - B*exp(-muz_cu)*exp(c_ba1))^CHI*((AA*exp(-vf_cup))/(exp(evf_cu)^(1/(ALFA - 1))*exp(muz_cup)^((CHI - 1)*(CHI0 - 1))))^ALFA)/(exp(c_cup) - B*exp(-muz_cup)*exp(c_cu))^CHI - 1;
f(29,1)=(BETTA*exp(-d_cu)*exp(-p20_cu)*exp(-pai_cup)*exp(d_cup)*exp(p19_cup)*exp(muz_cup)^(CHI*(CHI0 - 1) - CHI0)*(exp(c_cu) - B*exp(-muz_cu)*exp(c_ba1))^CHI*((AA*exp(-vf_cup))/(exp(evf_cu)^(1/(ALFA - 1))*exp(muz_cup)^((CHI - 1)*(CHI0 - 1))))^ALFA)/(exp(c_cup) - B*exp(-muz_cup)*exp(c_cu))^CHI - 1;
f(30,1)=(BETTA*exp(-d_cu)*exp(-p21_cu)*exp(-pai_cup)*exp(d_cup)*exp(p20_cup)*exp(muz_cup)^(CHI*(CHI0 - 1) - CHI0)*(exp(c_cu) - B*exp(-muz_cu)*exp(c_ba1))^CHI*((AA*exp(-vf_cup))/(exp(evf_cu)^(1/(ALFA - 1))*exp(muz_cup)^((CHI - 1)*(CHI0 - 1))))^ALFA)/(exp(c_cup) - B*exp(-muz_cup)*exp(c_cu))^CHI - 1;
f(31,1)=(BETTA*exp(-d_cu)*exp(-p22_cu)*exp(-pai_cup)*exp(d_cup)*exp(p21_cup)*exp(muz_cup)^(CHI*(CHI0 - 1) - CHI0)*(exp(c_cu) - B*exp(-muz_cu)*exp(c_ba1))^CHI*((AA*exp(-vf_cup))/(exp(evf_cu)^(1/(ALFA - 1))*exp(muz_cup)^((CHI - 1)*(CHI0 - 1))))^ALFA)/(exp(c_cup) - B*exp(-muz_cup)*exp(c_cu))^CHI - 1;
f(32,1)=(BETTA*exp(-d_cu)*exp(-p23_cu)*exp(-pai_cup)*exp(d_cup)*exp(p22_cup)*exp(muz_cup)^(CHI*(CHI0 - 1) - CHI0)*(exp(c_cu) - B*exp(-muz_cu)*exp(c_ba1))^CHI*((AA*exp(-vf_cup))/(exp(evf_cu)^(1/(ALFA - 1))*exp(muz_cup)^((CHI - 1)*(CHI0 - 1))))^ALFA)/(exp(c_cup) - B*exp(-muz_cup)*exp(c_cu))^CHI - 1;
f(33,1)=(BETTA*exp(-d_cu)*exp(-p24_cu)*exp(-pai_cup)*exp(d_cup)*exp(p23_cup)*exp(muz_cup)^(CHI*(CHI0 - 1) - CHI0)*(exp(c_cu) - B*exp(-muz_cu)*exp(c_ba1))^CHI*((AA*exp(-vf_cup))/(exp(evf_cu)^(1/(ALFA - 1))*exp(muz_cup)^((CHI - 1)*(CHI0 - 1))))^ALFA)/(exp(c_cup) - B*exp(-muz_cup)*exp(c_cu))^CHI - 1;
f(34,1)=(BETTA*exp(-d_cu)*exp(-p25_cu)*exp(-pai_cup)*exp(d_cup)*exp(p24_cup)*exp(muz_cup)^(CHI*(CHI0 - 1) - CHI0)*(exp(c_cu) - B*exp(-muz_cu)*exp(c_ba1))^CHI*((AA*exp(-vf_cup))/(exp(evf_cu)^(1/(ALFA - 1))*exp(muz_cup)^((CHI - 1)*(CHI0 - 1))))^ALFA)/(exp(c_cup) - B*exp(-muz_cup)*exp(c_cu))^CHI - 1;
f(35,1)=(BETTA*exp(-d_cu)*exp(-p26_cu)*exp(-pai_cup)*exp(d_cup)*exp(p25_cup)*exp(muz_cup)^(CHI*(CHI0 - 1) - CHI0)*(exp(c_cu) - B*exp(-muz_cu)*exp(c_ba1))^CHI*((AA*exp(-vf_cup))/(exp(evf_cu)^(1/(ALFA - 1))*exp(muz_cup)^((CHI - 1)*(CHI0 - 1))))^ALFA)/(exp(c_cup) - B*exp(-muz_cup)*exp(c_cu))^CHI - 1;
f(36,1)=(BETTA*exp(-d_cu)*exp(-p27_cu)*exp(-pai_cup)*exp(d_cup)*exp(p26_cup)*exp(muz_cup)^(CHI*(CHI0 - 1) - CHI0)*(exp(c_cu) - B*exp(-muz_cu)*exp(c_ba1))^CHI*((AA*exp(-vf_cup))/(exp(evf_cu)^(1/(ALFA - 1))*exp(muz_cup)^((CHI - 1)*(CHI0 - 1))))^ALFA)/(exp(c_cup) - B*exp(-muz_cup)*exp(c_cu))^CHI - 1;
f(37,1)=(BETTA*exp(-d_cu)*exp(-p28_cu)*exp(-pai_cup)*exp(d_cup)*exp(p27_cup)*exp(muz_cup)^(CHI*(CHI0 - 1) - CHI0)*(exp(c_cu) - B*exp(-muz_cu)*exp(c_ba1))^CHI*((AA*exp(-vf_cup))/(exp(evf_cu)^(1/(ALFA - 1))*exp(muz_cup)^((CHI - 1)*(CHI0 - 1))))^ALFA)/(exp(c_cup) - B*exp(-muz_cup)*exp(c_cu))^CHI - 1;
f(38,1)=(BETTA*exp(-d_cu)*exp(-p29_cu)*exp(-pai_cup)*exp(d_cup)*exp(p28_cup)*exp(muz_cup)^(CHI*(CHI0 - 1) - CHI0)*(exp(c_cu) - B*exp(-muz_cu)*exp(c_ba1))^CHI*((AA*exp(-vf_cup))/(exp(evf_cu)^(1/(ALFA - 1))*exp(muz_cup)^((CHI - 1)*(CHI0 - 1))))^ALFA)/(exp(c_cup) - B*exp(-muz_cup)*exp(c_cu))^CHI - 1;
f(39,1)=(BETTA*exp(-d_cu)*exp(-p30_cu)*exp(-pai_cup)*exp(d_cup)*exp(p29_cup)*exp(muz_cup)^(CHI*(CHI0 - 1) - CHI0)*(exp(c_cu) - B*exp(-muz_cu)*exp(c_ba1))^CHI*((AA*exp(-vf_cup))/(exp(evf_cu)^(1/(ALFA - 1))*exp(muz_cup)^((CHI - 1)*(CHI0 - 1))))^ALFA)/(exp(c_cup) - B*exp(-muz_cup)*exp(c_cu))^CHI - 1;
f(40,1)=(BETTA*exp(-d_cu)*exp(-p31_cu)*exp(-pai_cup)*exp(d_cup)*exp(p30_cup)*exp(muz_cup)^(CHI*(CHI0 - 1) - CHI0)*(exp(c_cu) - B*exp(-muz_cu)*exp(c_ba1))^CHI*((AA*exp(-vf_cup))/(exp(evf_cu)^(1/(ALFA - 1))*exp(muz_cup)^((CHI - 1)*(CHI0 - 1))))^ALFA)/(exp(c_cup) - B*exp(-muz_cup)*exp(c_cu))^CHI - 1;
f(41,1)=(BETTA*exp(-d_cu)*exp(-p32_cu)*exp(-pai_cup)*exp(d_cup)*exp(p31_cup)*exp(muz_cup)^(CHI*(CHI0 - 1) - CHI0)*(exp(c_cu) - B*exp(-muz_cu)*exp(c_ba1))^CHI*((AA*exp(-vf_cup))/(exp(evf_cu)^(1/(ALFA - 1))*exp(muz_cup)^((CHI - 1)*(CHI0 - 1))))^ALFA)/(exp(c_cup) - B*exp(-muz_cup)*exp(c_cu))^CHI - 1;
f(42,1)=(BETTA*exp(-d_cu)*exp(-p33_cu)*exp(-pai_cup)*exp(d_cup)*exp(p32_cup)*exp(muz_cup)^(CHI*(CHI0 - 1) - CHI0)*(exp(c_cu) - B*exp(-muz_cu)*exp(c_ba1))^CHI*((AA*exp(-vf_cup))/(exp(evf_cu)^(1/(ALFA - 1))*exp(muz_cup)^((CHI - 1)*(CHI0 - 1))))^ALFA)/(exp(c_cup) - B*exp(-muz_cup)*exp(c_cu))^CHI - 1;
f(43,1)=(BETTA*exp(-d_cu)*exp(-p34_cu)*exp(-pai_cup)*exp(d_cup)*exp(p33_cup)*exp(muz_cup)^(CHI*(CHI0 - 1) - CHI0)*(exp(c_cu) - B*exp(-muz_cu)*exp(c_ba1))^CHI*((AA*exp(-vf_cup))/(exp(evf_cu)^(1/(ALFA - 1))*exp(muz_cup)^((CHI - 1)*(CHI0 - 1))))^ALFA)/(exp(c_cup) - B*exp(-muz_cup)*exp(c_cu))^CHI - 1;
f(44,1)=(BETTA*exp(-d_cu)*exp(-p35_cu)*exp(-pai_cup)*exp(d_cup)*exp(p34_cup)*exp(muz_cup)^(CHI*(CHI0 - 1) - CHI0)*(exp(c_cu) - B*exp(-muz_cu)*exp(c_ba1))^CHI*((AA*exp(-vf_cup))/(exp(evf_cu)^(1/(ALFA - 1))*exp(muz_cup)^((CHI - 1)*(CHI0 - 1))))^ALFA)/(exp(c_cup) - B*exp(-muz_cup)*exp(c_cu))^CHI - 1;
f(45,1)=(BETTA*exp(-d_cu)*exp(-p36_cu)*exp(-pai_cup)*exp(d_cup)*exp(p35_cup)*exp(muz_cup)^(CHI*(CHI0 - 1) - CHI0)*(exp(c_cu) - B*exp(-muz_cu)*exp(c_ba1))^CHI*((AA*exp(-vf_cup))/(exp(evf_cu)^(1/(ALFA - 1))*exp(muz_cup)^((CHI - 1)*(CHI0 - 1))))^ALFA)/(exp(c_cup) - B*exp(-muz_cup)*exp(c_cu))^CHI - 1;
f(46,1)=(BETTA*exp(-d_cu)*exp(-p37_cu)*exp(-pai_cup)*exp(d_cup)*exp(p36_cup)*exp(muz_cup)^(CHI*(CHI0 - 1) - CHI0)*(exp(c_cu) - B*exp(-muz_cu)*exp(c_ba1))^CHI*((AA*exp(-vf_cup))/(exp(evf_cu)^(1/(ALFA - 1))*exp(muz_cup)^((CHI - 1)*(CHI0 - 1))))^ALFA)/(exp(c_cup) - B*exp(-muz_cup)*exp(c_cu))^CHI - 1;
f(47,1)=(BETTA*exp(-d_cu)*exp(-p38_cu)*exp(-pai_cup)*exp(d_cup)*exp(p37_cup)*exp(muz_cup)^(CHI*(CHI0 - 1) - CHI0)*(exp(c_cu) - B*exp(-muz_cu)*exp(c_ba1))^CHI*((AA*exp(-vf_cup))/(exp(evf_cu)^(1/(ALFA - 1))*exp(muz_cup)^((CHI - 1)*(CHI0 - 1))))^ALFA)/(exp(c_cup) - B*exp(-muz_cup)*exp(c_cu))^CHI - 1;
f(48,1)=(BETTA*exp(-d_cu)*exp(-p39_cu)*exp(-pai_cup)*exp(d_cup)*exp(p38_cup)*exp(muz_cup)^(CHI*(CHI0 - 1) - CHI0)*(exp(c_cu) - B*exp(-muz_cu)*exp(c_ba1))^CHI*((AA*exp(-vf_cup))/(exp(evf_cu)^(1/(ALFA - 1))*exp(muz_cup)^((CHI - 1)*(CHI0 - 1))))^ALFA)/(exp(c_cup) - B*exp(-muz_cup)*exp(c_cu))^CHI - 1;
f(49,1)=(BETTA*exp(-d_cu)*exp(-p40_cu)*exp(-pai_cup)*exp(d_cup)*exp(p39_cup)*exp(muz_cup)^(CHI*(CHI0 - 1) - CHI0)*(exp(c_cu) - B*exp(-muz_cu)*exp(c_ba1))^CHI*((AA*exp(-vf_cup))/(exp(evf_cu)^(1/(ALFA - 1))*exp(muz_cup)^((CHI - 1)*(CHI0 - 1))))^ALFA)/(exp(c_cup) - B*exp(-muz_cup)*exp(c_cu))^CHI - 1;
f(50,1)=exp(-c_ba1p)*exp(c_cu) - 1;
f(51,1)=RHOz*log(exp(muz_cu)/MUZss) - log(exp(muz_cup)/MUZss);
f(52,1)=RHOd*log(exp(d_cu)) - log(exp(d_cup));
f(53,1)=RHOn*log(exp(n_cu)) - log(exp(n_cup));
f(54,1)=RHOp*log(exp(paistar_cu)) - log(exp(paistar_cup));
f(55,1)=RHOa*log(exp(a_cu)) - log(exp(a_cup));
% END DISPLAYING f
        Euler(:,1) = Euler(:,1) + wGrid(i,1)*f(:,1); 
    end 
    ResEuler(k,:) = transpose(Euler(:,1));
end 
delete(gcp('nocreate'))